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Abstract 

We use the Litvinov-Maslov correspondence principle to reduce and hybridize networks of biochemical 
reactions. We apply this method to a cell cycle oscillator model. The reduced and hybridized model can be 
used as a hybrid model for the cell cycle. We also propose a practical recipe for detecting quasi-equilibrium 
QE reactions and quasi-steady state QSS species in biochemical models with rational rate functions and 
use this recipe for model reduction. Interestingly, the QE/QSS invariant manifold of the smooth model 
and the reduced dynamics along this manifold can be put into correspondence to the tropical variety of 
the hybridization and to sliding modes along this variety, respectively 
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1 Introduction. 

Systems biology develops biochemical dynamic models of various cellular processes such as signalling, 
metabolism, gene regulation. These models can reproduce complex spatial and temporal dynamic be- 
havior observed in molecular biology experiments. In spite of their complex behavior, currently available 
dynamical models are relatively small size abstractions, containing only tens of variables. This modest 
size results from the lack of precise information on kinetic parameters of the biochemical reactions on one 
hand, and of the limitations of parameter identification methods on the other hand. Further limitations 
can result from the combinatorial explosion of interactions among molecules with multiple modifications 
and interaction sites |DFF + 07] . In middle out modeling strategies small models can be justified by saying 
that one looks for an optimal level of complexity that captures the salient features of the phenomenon 
under study. The ability to choose the relevant details and to omit the less important ones is part of 
the art of the modeler. Beyond modeler's art, the success of simple models relies on an important prop- 
erty of large dynamical systems. The dynamics of multiscale, dissipative, large biochemical models, can 
be reduced to that of simpler models, that were called dominant subsystems (RGZL08 I IGRZlO l fGR08 . 
Simplified, dominant subsystems contain less parameters and are more easy to analyze. The choice of 
the dominant subsystem depends on the comparison among the time scales of the large model. Among 
the conditions leading to dominance and allowing to generate reduced models, the most important are 
quasi-equilibrium (QE) and the quasi-steady state (QSS) approximations [GRZlfJ . In nonlinear systems, 
timescales and together with them dominant subsystems can change during the dynamics and undergo 
more or less sharp transitions. The existence of these transitions suggests that a hybrid, discrete/continous 



framework is well adapted for the description of the dynamics of large nonlinear systems with multiple 
time scales |UDR09l INVRlOl INVRllj . 

The notion of dominance can be exploited to obtain simpler models from larger models with multiple 
separated timescales and to assemble these simpler models into hybrid models. This notion is asymptotic 
and a natural mathematical framework to capture multiple asymptotic relations is the tropical geometry. 
Motivated by applications in mathematical physics [LM96 , systems of polynomial equations [Stu02], etc., 
tropical geometry uses a change of scale to transform nonlinear systems into discontinuous piecewise linear 
systems. The tropicalization is a robust property of the system, remaining constant for large domains of 
parameter values; it can reveal qualitative stable features of the system's dynamics, such as various types 
of attractors. Thus, the use of tropicalization to model large systems in molecular biology could be a 
promising solution to the problem of incomplete or imprecise information on the kinetic parameters. 

In this paper we propose a method for reduction and hybridization of biochemical networks. This 
method, based on tropical geometry, could be used to automatically produce the simple models that are 
needed in middle-out approaches of systems biology. 

2 Biochemical networks with rational rate functions. 

Systems biology models use the formalism of chemical kinetics to model dynamics of cellular processes. 
We consider here that the molecules of various species are present in sufficient large numbers and that 
stochastic fluctuations are negligible as a consequence of the law of large numbers and/or of the averag- 
ing theorem [CDR09] . We also consider that space transport phenomena are sufficiently rapid such that 
the well stirred reactor hypothesis is valid. In these conditions, the dynamics of the biochemical system 
can be described by systems of differential equations. In chemical kinetics, enzymatic reactions are of- 
ten presented as indivisible entities characterized by stoichiometry vectors and rate functions. However, 
each enzymatic reaction can be decomposed into several steps that define the reaction mechanism. The 
resulting stoichiometry and global rate depend on the mechanism. Several methods were designed for 
calculating effective rates of arbitrarily complex mechanisms. For linear mechanisms King and Altman 
[KA56| proposed a graphical method to compute global rates; these are rational functions of the con- 
centrations (an example is the Michaelis-Menten equation). Yablovsky and Lazman |LY08bj studied the 
same problem for non-linear mechanisms and found that in this case the reaction rates are solutions of 
polynomial equations; these can be solved by radicals in a few number of cases and can be calculated by 
multi-variate hypergeometric series in general [LY08b . Truncation of these series to finite order leads to 
rational approximations of the reaction rates. 

In chemical kinetics with rational reaction rates the concentration Xi of the i-th component follows the 
ordinary differential equation: 

^=P i (x)/Q i (x), (1) 

where Pi(x) = X^q<ea a i.aX a , Qi{x) = X^es bi^x @ , are polynomials and we have 1 < i < n. Here 
x a — x" 1 x% 2 . . . x™ n , x@ = x^ 1 x^ 2 ■ ■ ■ x^ n , ai ia ,bi y p, are nonzero real numbers, and A^Bi are finite 
subsets of N™ called supports of Pi and Qi. 

A simple example of model with rational reaction rates is the minimal cell cycle oscillator model 
proposed by Tyson |Tys91| . This example will be studied throughout the paper. The dynamics of this 
nonlinear model that contains 5 species and 7 reactions is described by a system of 5 polynomial differential 
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equations: 

y'x = k 9 y 2 - k 8 Vi + hy 3 , 

y' 2 = hyi - k 9 y 2 - k 3 y 2 y 5 , 

V 3 = k' 4 y 4 + k 4 y 4 y%/C 2 ~ k e y 3 , 

y'i = -k' 4 y 4 ~ k 4 y 4 yj/C 2 + k 3 y 2 y 5 , 

2/5 =ki- k 3 y 2 y 5 , (2) 

where y\ + y 2 + y 3 + y 4 = C. 



3 Hybridization and tropical geometry. 

Tropical geometry is a new branch of algebraic geometry that studies the asymptotic properties of vari- 
eties. While algebraic geometry deals with polynomial functions, tropical geometry deals with piecewise 
linear functions with integer directing slopes. Tropical geometry has a growing number of applications in 
enumerative problems in nonlinear equation solving R oj03| , statistics [P"S04 , traffic optimization |AublO| . 

The logarithmic transformation Ui = logxi, 1 < i < n, well known for drawing graphs on logarithmic 
paper, plays a central role in tropical geometry |Vir08] . By abus de langage, here we call logarithmic paper 
the image of K™ by the logarithmic transformation, even if n > 2. Monomials M(x) — a a x a with positive 
coefficients a a > 0, become linear functions, logM — loga a + < a,log(x) >, by this transformation. 
Furthermore, the euclidian distance on the logarithmic paper is a good measure of separation (see next 
section). 

Litvinov and Maslov [LMS011 ILM96j proposed a heuristic (correspondence principle) allowing to trans- 
form mathematical objects (integrals, polynomials) into their quantified (tropical) versions. According to 
this heuristic, to a polynomial with positive real coefficients X^ae^t a a& a > one associates the max-plus 
polynomial max a ^A{log(a a ) + < log(x),a >}. 

We adapt this heuristic to associate a piecewise-smooth hybrid model to the systems of rational ODEs 

0. 

Definition 1 

We call tropicalization of the smooth ODE system the following piecewise-smooth system: 
dx ' 

-— ^ = Siexp[max aeAi {log(\ai, a \)+ <u,a>} - max l 3 eB A lo g(\ b i,p\)+ <u,(3 >}], (3) 

where u = {logxx, . . . ,logx n ), s 4 = sign(a itamax )sign(b hl 3 ma;i: ) and a 4 , Qmax , a max € A, (respectively, 
bitPmaa, > Pmax € Bi) denotes the coefficient of a monomial of the numerator (respectively, of the denomi- 
nator) for which the maximum occurring in ^ is attained. 

In a different notation this reads: 

-J^ = Dom{ai^x a } aeA jDom{bi,pX P } aeBz , (4) 

where Dom{ai^ a x a } ae Ai = sign(ai lCtmaa .)exp[rnax a ^A i {log(\ai ta \)+ <u,a >}]. 
Finally, the tropicalization can be written with Heaviside functions: 

dxj _ EagA, a i,aX a Y\ a ,^ a e{< a - a',log(x) > +log(\ai, a \) - log{\ai, a >\)) 
dt ~ E^ eSl b itP xP l\p>^ 0(< P - /?', log(x) > +log(\b l . f} \) - log(\b hP ,\)) ' 

where 9(x) = 1 if x > 0, if not. 

The following definitions are standard and will be used throughout the paper: 
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Definition 2 

The Newton polytope of a polynomial P(x) — J2 a <=A a aX a J 's defined as the convex hull of the support of 
P, New(P) = conv{A). 

Definition 3 

The max-plus polynomial P T (x) — max{log\a a \+ < a,log(x) >} is called the tropicalization of P(x). 
The logarithmic function is dehned as log(x) : K™ — > E", log(x)i = log(xi). 

Definition 4 

The set of points x £ K™ where P T (x) is not smooth is called tropical variety. Alternative names are used 
such as logarithmic limit sets, Bergman fans, Bieri-Groves sets, or non Archimedean amoebas \PT05l 

In two dimensions, a tropical variety is a tropical curve made of several half-lines (tentacles) and finite 
intervals [Mikj . A tropical line corresponds to only three monomials and is made of three half lines sharing 
a common point. The tentacles and the intervals of the tropical variety are orthogonal to the edges and 
point to the interior of the Newton polygon |PT05| (see Figjl]). 

4 Dominance and separation. 

The above heuristic is related to the notion of dominance. Actually we have replaced each polynomial in 
the rational function by the dominant monomial. Dominance of monomials has an asymptotic meaning 
inside cones of the logarithmic paper. For instance x a dominates x" on the half plane < log(x),a—/3 >> 
of the logarithmic paper. We have x$ jx a — > when the limit is taken along lines in this half plane. 

For practical applications, we would also need a finite scale notion of dominance. 

Let Mi{x) = a ai x ai and M^x) = a a2 x a2 be two monomials. We define the following binary relations: 

Definition 5 (Separation) 

Mi and are separated on a domain D C i?" at a level p > if \log(\a ai \x ai ) — log(\a a2 |a; Q2 )| > p for 
all x e D. 

On logarithmic paper, two monomials are separated on the domain D, if D is separated by the euclidian 
distance p from the hyperplane < log(x), ol\ — o>2 >— log\a a2 \ — log\a ai |. 

Definition 6 (Dominance) 

The monomial M\ dominates the monomial at the level p > 0, Mj > p M2, if log(\a ai \x ai ) > 
log(\a a2 \x a2 ) + p for all x e D C E™. 

Dominance is a partial order relation on the set of multivariate monomials defined on subsets of M™ . 

5 Dominance and global reduction of large models. 

There are two simple methods for model reduction of nonlinear models with multiple timescales: the 
quasi-equilibrium (QE) and the quasi-steady state (QSS) approximations. As discussed in [GRZ10J, these 
two approximations are physically and dynamically distinct. Here we present a method allowing to detect 
QE reactions and QSS species. 

Like in [RGZL08 , the first step of the method is to detect the "slaved" species, i.e. the species that 
obey quasi-steady state equations. These can be formally defined by introducing the notion of imposed 
trace. Given the traces x(t) of all the species, the imposed trace of the i-th species is a real solution x*(t) 
of the polynomial equation Pi(xi(t), . . . , Xi^\{t), x*(t), Xi + i(t), . . . , x n {t)) = 0. Eventually, there may be 
several imposed traced, because a polynomial equation can have several real solutions. 
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Definition 7 

We say that a species is slaved if the distance between the traces Xi (t) and some imposed trace x* (t) is 
small on some interval, suptei\log(xi(t)) — log(x*(t))\ < S, for some S > sufficiently small. A species is 
globally slaved if I = (T, oo) for some T > 0. 

Slaved species are good candidates for QSS species and this criterion was used to identify QSS species 
in |RGZL08] . More generally, slaved species are involved in rapid processes, but are not always QSS. 
Actually, two distinct cases lead to slaved species. 

Quasi- equilibrium. A system with fast, quasi-equilibrium reactions has the following structure GRZ10 : 



= J2 R^h s + - E RA*h f , (6) 



dx 

s,slow j, fast 



where e > is a small parameter 7 s , 7^ e Z™ are stoichiometric vectors. The reaction rates R s (x), Rf(x) 
are considered rational functions of x. 

To separate slow/fast variables, we have to study the spaces of linear conservation law of the initial 
system (|6| and of the following fast subsystem: 

fjast 

In general, the system (JfiJ can have several conservation laws. These are linear functions b 1 (x) , . . . , b m (x) 
of the concentrations that are constant in time. The conservation laws of the system ([T]) provide vari- 
ables that are constant on the fast timescale. If they are also conserved by the full dynamics, the system 
has no slow variables (variables are either fast or constant). In this case, the dynamics of the fast 
variables is simply given by Eq.([7|. Suppose now that the system ^ has some more conservation laws 
b m+1 (x), . . . , b m+l (x), that are not conserved by the full system Then, these provide the slow variables 
of the system. The fast variables are those x$ such that (7* 7^ 0, for some fast reaction /. 

Let us suppose that the fast system ^ has a stable steady state that is a solution of the QE equations 
(augmented by the conservation laws of the fast system): 

Z f ,fastRA x h f = °> ( 8 ) 

b l {x)^C ll l<i<m + l. (9) 

By classical singular perturbation methods |TVS 85 , Was 65] one can show that the fast variables can 
be decomposed as Xi = £i + r\i where Xi satisfy the QE equations pi and 77, = 0(e), meaning that the 
fast variables slaved iGRZlOj . 

Let Pi, Pi be the numerators of the rational functions J^s slow Rs( x )~fi + (E/ fast Rf( x )li an d 
J2f fast Rf( x )'~fi' respectively. We call Pj the pruned version of Pj. When e is small enough, the monomials 
of the pruned version Pj dominate the monomials of Pj. This suggests a practical recipe for identifying 
QE reactions: 
Algorithm 1 

Step 1: Detect slaved species. 

Step 2: For each Pj corresponding to slaved species, compute the pruned version Pj by eliminating all 
monomials that are dominated by other monomials of Pj. 

Step 3: Identify, in the structure of Pj the forward and reverse rates of QE reactions. This step could be 
performed by recipes presented in [SHlOj . 



■5 



Quasi-steady state. In the most usual version of QSS approximation SS89J, the species are split in 
two groups with concentration vectors x s ("slow" or basic components) and x^ ("fast" or QSS species). 

Quasi-steady species (also called radicals or fast intermediates) are low-concentration, slaved species. 
Typically, QSS species are consumed (rather than produced) by fast reactions. The small parameter e 
used in singular perturbation theory is now the ratio of small concentrations of fast intermediates to the 
concentration of other species. After rescaling x s and x' to order one, the set of kinetic equations reads: 



dx s 
dx* 



= W s (x s ,xf), (10) 
= (l/e)Wf(x s ,xf), (11) 
where the functions W s , W* and their derivatives are of order one (0 < e << 1). 



Let us suppose that the fast dynamics (11) has a stable steady state. The standard singular pertur- 
bation theory [TVS851 IWas65| provides the QSS algebraic condition W^(x s , x$) = which means that 
fast species x^ are slaved. These equations, together with additional balances for x? (conservation laws) 
are enough to deduce the fast variables x? as functions of the slow variables x s and to eliminate them 



[YBGE91I lLY0"8al IRGZL08] . The slow dynamics is given by Eq.([T0|. 

In networks with rational reaction rates the components of W* (x s ,x?) are rational functions. Like 
for QE we can define Pi as numerators of Wf . The difference between QSS conditions with respect to 
QE situation is that in the pruned polynomial Pi one can no longer find forward and backward rates of 
QE reactions, ie the step 3 of Algorithm 1 will not identify reversible reactions. Alternatively, one can 
realize that slaved species can have relatively large concentrations, in which case they are not QSS species. 
However, it is difficult to say which concentration value separates QSS from non QSS species among slaved 
species, hence the former, qualitative criterion is better. 



6 Sliding modes of the tropicalization. 

A notable phenomenon resulting from tropicalization is the occurrence of sliding modes. Sliding modes 
are well known for ordinary differential equations with discontinuous vector fields FA88 . In such systems, 
the dynamics can follow discontinuity hypersurfaces where the vector field is not defined. 

The conditions for the existence of sliding modes are generally intricate. However, when the discon- 
tinuity hypersurfaces are smooth and n — 1 dimensional (n is the dimension of the vector field) then the 
conditions for sliding modes read: 

< n + (x), f + (x) >< 0, < ri-(x),f-(x) >< 0, x G S, (12) 

where /+, /_ are the vector fields on the two sides of S and n + = — n_ are the interior normals. 

Let us consider that the smooth system ([!]) has quasi-steady state species or quasi-equilibrium reactions. 
In this case, the fast dynamics reads: 

dr- 1 ~ 

= -Pi(x)/Qi(x), i fast, (13) 

where Pi(x), Qi(x) are pruned versions of Pi, Qi, and e is the small, singular perturbation parameter. 
For sufficiently large times, the fast variables satisfy (to O(e)): 

Pi(x)=0, i fast. (14) 

The pruned polynomial is usually a fewnomial (contains a small number of monomials). In particular, 
let us consider the case when only two monomials remain after pruning, Pi(x) = a\X ai +a2X a2 . Then, the 
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equation ( 14) defines a hyperplane S — {< log(x), a\~a 2 >— Zo.g( |ai | / 1 a 2 1 )}. This hyperplane belongs to 
the tropical variety of Pi, because it is the place where the monomial x ai switches to x a2 in the max-plus 
polynomial defined by Pj. For e small, the QE of QSS conditions guarantee the existence of an invariant 
manifold A4 e , whose distance to S is 0(e). 

Let n + ,n_ defined as above and let (f + ) l = -^- ) a 1 x ai , (/_), = -^^a 2 x a2 , f t = -[(/+)* + (/_)*] for 

i fast, (f+)j = (f-)j = fj = 7t , for j not fast. Then, the stability conditions for the invariant manifold 

Qj 

read < n + (x + ), f(x + ) >< 0, < n_(.T_), f(x-)) >< 0, where x +1 x^ are close to M e on the side towards 
which points n + and n_, respectively. We note that \{f+)i(x + )\ > \(f-)i(x + )\. Thus, < n + , f >= 

l A n +)i[(f+h + + E J ,notfaJn+Uf+)i and < n+,/ + >= i(n + ),(/ + ), + Ei,nat/„t(n+)i(/+)i- 

Thus, if < n + , / >< 0, then for e small enough (n + )i(/ + )j < and < n + ,f + >< because < n + , / >>< 
n + , / + >. Similarly, we show that < n_, / >< implies < ra_, /_ >< 0. This proves the following 

Theorem 1 

If the smooth dynamics obeys QE or QSS conditions and if the pruned polynomial P defining the fast 
dynamics is a 2-nomial, then the QE or QSS equations define a hyperplane of the tropical variety of P. 
The stability of the QE of QSS manifold implies the existence of a sliding mode of the tropicalization 
along this hyperplane. 

The converse result, i.e. deducing the stability of the QE/QSS manifold from the existence of a sliding 
mode on the tropical variety may be wrong. Indeed, it is possible for a trajectory of the smooth system to 
be close to a hyperplane of the tropical variety carrying a sliding mode and where the QE/QSS equations 
are satisfied identically. However, as we will see in the next section, this trajectory can leave the hyperplane 
sooner than the sliding mode. 



7 From smooth to hybrid models via reduction. 

Starting with the system ^ we first reduce it to a simpler model. The analysis of the model is performed 
for the values of parameters from |Tys91| , namely k\ = 0.015, = 200, = 180, k' 4 = 0.018, k§ = 1, = 
0.6, k$ = 1000000, fc 9 = 1000; 

In order to do that we generate one or several traces (trajectories) Di(t). The smooth system has a 
stable periodic trace which is a limit cycle attractor. We also compute the imposed traces y* (t) that are 
solutions of the equations: 

&93/2(*) - hy*(t) + k & y z {t) = 0, 
k 8 yi(t) - kgylif) - k 3 y 2 (t)y 5 (t) = 0, 
Ky 4 (t) + k m {t)y* 3 2 (t)/C 2 - k e y* 3 {t) = 0, 
-KvW) ~ hyl(t)y 2 3 (t)/C 2 + k 3 y 2 (t)y 5 (t) = 0, 

ki - k 3 y 2 (t)y* 5 (t) = 0. (15) 

We find that, for three species 2/1,1/2, and j/5, the distance between the traces y*(t) and yi(t) is small 
for all times which means that these species are slaved on the whole limit cycle (Figure [T] top left) . Also, 
we have a global conservation law yi + y 2 + y 3 + 2/4 = C, that can be obtained by summing the first four 
differential equations in ([2]). The three quasi-steady state equations for the three slaved species have to 
be solved jointly with the global conservation law: 

fcgy2 - k 8 yi + k 6 y 3 = 0, 
ksVi - k 9 y 2 - k 3 y 2 y 5 = 0, 
ki - k 3 y 2 y 5 = 0, 
Vi + 2/2 + 2/3 + 2/4 = C. (16) 
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Figure 1: (top left) Detection of slaved species by comparing traces to imposed traces: the species yi, y%, y§ 
are slaved globally, the species 3/3,2/4 are slaved on intervals Qz,Qi, respectively, (top right) Comparison 
of monomials of the polynomial systems of quasi-steady state equations, (bottom left) Newton polygons 
and inner normals of the reduced two dimensional polynomial model, (bottom right) Phase portrait on 
logarithmic paper of the reduced two dimensional model. We represent the two tropical curves (the tripods 
graphs, a red and a blue one), the modes (smooth vector fields within domains bordered by tropical curves 
tentacles), the smooth and tropicalized limit cycles. The tropicalized cycle contains two sliding modes 
53,6*4 corresponding to the intervals Q3, Q4 on which y 3 , y 4 are quasi-stationary, respectively. 



Comparison of the monomials (for values of parameters as above) in this system shows that max(k%y\, fcg?^) >~ 
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keys, and max(kg,yx, fegj/2) >- k 3 y 2 y§ (FigjTjtop right) which leads to the pruned system: 

hyi ~ k 9 V2 = 0, 
hyi ~ hy 2 = 0, 
ki - k 3 y 2 y 5 = 0, 
2/i + 2/2 + 2/3 + 2/4 = C. (17) 

The first two equations are identical and correspond to quasi-equilibrium of the reaction between yi 
and y 2 . The third equation means that 7/5 is a quasi-steady state species. The pruned system allows the 
elimination of the variables 2/1:2/2,2/5- The slow variable 2/12 =2/1+2/2 demanded by the quasi-equilibrium 
condition (this is a conservation law of the fast system) can be eliminated by using the global conservation 
law. 

We note that the dominance relations leading to the pruned equations were found numerically in a 
neighborhood of the periodic trace. This means that QE and QSS approximations are valid at least on the 



limit cycle. More global testing of these relations will be presented elsewhere. Note that the system ( 16 ) 



can be solved also without pruning. However, ( 16 1 has four independent equations allowing to eliminate 



four of the five dynamic variables leading to a one dimensional dynamical system. It turns out that the 



correct application of the QE and QSS approximations has to use (17 1 and not (16 1. 

After elimination, we obtain the following reduced differential-algebraic dynamical system: 

2/3 = k' 4 y 4 + k 4 y 4 y%/C 2 - k 6 y 3 , 
2/4 = -k' 4 y 4 - k 4 y 4 yl/C' 2 + hi, 
2/1 = (C- 2/3 -2/4)^9/(^8 + ^9), 
2/2 = (C - 2/3 - 2/4)W(fcs + ko,), 

2/5 = h(k 8 +k 9 )/(k 3 k 8 (C -y 3 -y 4 ). (18) 
Now we tropicalize this reduced system. The tropicalization could have been done on the initial system 



in which case the pruned equations ( 17) would indicate that the reduced dynamics is a sliding mode of the 
tropicalized system on the two dimensional hypcrsurface fegj/i = k 9 y 2 , fei = k 3 y 2 y 5 , 2/1 + 2/2 + 2/3 + 2/4 = C . 
However, although the result (concerning the dynamics on the QE/QSS manifold) should be the same, it 



is much handier to tropicalize the reduced system ( 18 ). Indeed, the tropicalization of the full 5D system is 



difficult to visualize and would also produce complex modes that can not be reduced to 2D (these modes 
describe the fast relaxation to the QE/QSS manifold). 
The resulting hybrid model reads: 



(19) 



y' 3 = Dom{k 4 y 4 , k 4 y 4 y 2 /C 2 , -k 6 y 3 }, 
y 4 = Dom{-k' 4 y 4 , -k 4 y 4 y 2 /C 2 , ki}, 

or equivalently using Heaviside functions: 

l/ 3 = Kv^i-hi - 2u 3 )9(h 2 +u 4 - u 3 ) + ^ViVl6{hi + 2u 3 )d{h 1 +h 2 +u 4 + u 3 ) 

-k 6 y 3 0(-h 2 -u 4 + u 3 )0(-hi - h 2 -114- u 3 ), 
y 4 = -k' 4 y 4 6(-h 3 - 2u 3 )9{-h 4 + u 4 ) - ^y±ylQ{h 3 + 2u 3 )9(h 3 -h 4 + 2u 3 + u 4 ) 

ki9(h 4 - u 4 )6(-h 3 + h 4 - 2u 3 - u 4 ), (20) 

where hi = h 3 = log(k 4 /(k' 4 C 2 )), h 2 — log(k' 4 /k 6 ), h 4 = log(ki/k' 4 ). 

The Newton polygons of the polynomials k' 4 y 4 + k 4 y 4 y 2 /C 2 — k§y 3 and —k' 4 y 4 — k 4 y 4 y 2 /C 2 — k§y 3 are 
triangles (Figfl] bottom left). The two triangles share a common edge which is a consequence of the fact 
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that the reduced model have two reactions each one acting on the two species. The tentacles of the two 
tropical curves (in red and blue in Fig[l] bottom right) point in the same directions as the inner normals 
to the edges of the Newton polygons (the corresponding equations are hi + 2u^ — 0, hi + u 4 — u 3 = 0, 
h\ + hi + u 4 + M3 = for one and h 3 + 2u 3 = 0, /14 + u 4 = 0, h 3 — /14 + 2^3 + U4 = for the other). These 
tentacles (half lines) decompose the positive quarter plane into 6 sectors corresponding to the 6 modes of 
the hybrid model. 

In Fig [l] bottom right we have also represented the phase portrait of the reduced model on logarithmic 
paper. The dynamical variables are U3 = log{y 3 ) and 1*4 = log(yi). The vector field corresponding to 
u 3 = 2/3/2/3 an d u' 4 = 2/4/2/4 was computed with the dominant monomials in each plane sector as follows: 

u 4 — — kiys, u' 3 ~ ~ ^6 for the mode 1, 

u 4 = — k/iy 3 ,u 3 = ^42/32/4 for the mode 2, 

u'4 — kiyl ,u' 3 = k^yzyi for the mode 3, 

u 4 = kiy~[ ,u' 3 = k'^y4,y 3 for the mode 4, 
u 4 — k\y~l , u' 3 — —k% for the mode 5, 

u 4 = — k'^,v! 3 = k 4 y4y 3 1 for the mode 6. (21) 

Like the smooth system, the tropicalization has a stable periodic trajectory (limit cycle). This is 
represented together with the limit cycle trajectory of the smooth system in Fig[l] bottom right. The 
period of the tropicalized limit cycle is slightly changed with respect to the period of the smooth cycle. 
However, we can modulate the period of the tropicalized cycle and make it fit the period of the smooth 
cycle by acting on the moments of the mode change. This stands to displacing the tentacles of the tropical 
varieties parallel to the initial positions or equivalently, to changing the parameters hi, hi, h 3 , /14 while 
keeping hi = /13 which is a symmetry of the problem. 

The tropicalized system has piecewise smooth hybrid dynamics. Typically, it passes from one type 
of smooth dynamics (mode) described by one set of differential equations to another smooth dynamics 



(mode) described by another set of differential equations (the possible modes are listed in Eq.(21|). The 



command to change the mode is intrinsic and happens when the trajectory attains the tropical curve. 



However, if the sliding mode condition (12 1 is fulfilled the trajectory continues along some tropical curve 
tentacle instead of changing plane sector and evolve according to one of the interior modes ( plj ). The 
tropicalized limit cycle has two sliding modes (S4 and 53 in Fig[T]). The first one is along the half-line 
h 3 — h 4 + 2u 3 + u 4 = on the logarithmic paper (tentacle S4 on the red tropical curve in Fig{l]) . In order to 
check (12) we note that / + = {kiy± , — ko), f~ = (— fejj/f, —k e ), n + = —rT = (—1, —2). We have a sliding 



mode if —kiy^ 1 + 2ke < 0, meaning that the exit from the sliding mode occurs when M4 > log(ki/(2ke)). 
The second sliding mode is along the tentacle hi + U4 — U3 = (S 3 on the blue tropical curve in Fig[I]) . We 
have / + = (kiy~[ ,—ke), f~ = (kiy 4 1 ,k' 4 y 4 y 3 1 ), n + = —rT = (—1,1). The conditions (12) are fulfilled 



when kiy 4 x — k' 4 y 4 y 3 ~ 1 < which is satisfied on the entire tentacle. The exit from this second mode occurs 
at the end of the blue tropical curve tentacle. Interestingly, the sliding modes of the tropicalization can 
be put into correspondence with places on the smooth limit cycle where the smooth limit cycle acquires 
new QSS species. This can be seen in Figjljtop left. The species 2/3 becomes quasi-stationary on time 
intervals Q 3 that satisfy (with good approximation) the relation hi + u 4 — u 3 = and correspond to the 
sliding mode on the blue tropical curve. Also, the species 2/4 becomes quasi-stationary on very short time 
intervals Q4 that satisfy h 3 — h 4 + 2u 3 + u 4 = and correspond to the sliding mode on the red tropical 
curve. As pointed out in the preceding section, the trajectories of the smooth dynamics can evolve close 
to the tentacles, but leave them sooner than the sliding modes. 

We end this section with a study of the bifurcations of the ODE model and of its tropicalization. It is 
easy to check that there is only one degree of freedom describing the relative position of the two tropical 
curves. This is the distance between the origins of the tropical curves, that is given by the combination 
kik 4 k 4 kg . Thus, by changing any one of the parameters ki,k' 4 ,k4,ke we can invert the relative 



10 



position of the tropical curves and change the partition of the logarithmic paper into domains. This leads 
to two Hopf bifurcations of the ODE model and also two Hopf bifurcations of the tropicalization. The 
bifurcation of the tropicalization is discontinuous and can also be delayed with respect to the continuous 
bifurcation of the ODE model (Fig. 2). 




Figure 2: Hopf bifurcations of the smooth and tropicalized system, (left) The relative positions of the 

/ 1/2 1/2 1 

tropical curves can be changed by changing the combination k\k A k 4 k e . The first Hopf bifurcation 

/ 1/2 1/2 i 

corresponds to fcifc 4 k 4 k^ = 1, i.e. log(ki) — —4.61, when the tropical curves intersect in a single 
point. For the second Hopf bifurcation the relative position of the two tropical curves is no longer 
exceptional; the position of the bifurcation results from sliding modes stability analysis, (right) Amplitudes 
of oscillation are shown for the tropicalization (red) and for the smooth system (blue); 



8 Solving ordinary differential equations in triangular form 

We give a digest of a general algorithm for solving systems of the type ([T]) and more generally, an arbitrary 
system of ordinary differential equations: 

G j (x 1 ,x?\...,x[,x 2 ,x£\...,x£\...,x n ,x£\...,x£\t) = 0, 1 < j < N, (22) 

where Gj are differential polynomials of the order at most r in the derivatives x\ = d s Xi/dt s , s < r. 
Let the degrees of the differential polynomials Gj do not exceed d. Finally, for algorithmic complexity 
purposes we assume that the coefficients of Gj are integers with absolute values less than 2 l , the latter 
means that the bit-size of the coefficients l(Gj) < I 



In [Sei56 an algorithm was designed which works not only for ordinary differential systems like ( 22 ) , but 



even for systems of partial differential equations. For ordinary systems ( 22 ) the algorithm was improved in 
|Gri89| , although still its complexity is rather big (see below) . We describe the ingredients of the output 
(which has a triangular form) of the latter improved algorithm and provide for it the complexity bounds. 

The algorithm executes the consecutive elimination of the indeterminates x n , . . . , x\. The algorithm 
yields a partition P — {Pi}i<i<M of the space of the possible functions x\. Each Pi is given by a system of 
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an equation f% t \{x-\_,t) = and an inequality gi^(x\,t) 7^ for suitable differential polynomials /^i, g^\. 
Then the algorithm yields an equation fi t 2(xi,xa,t) = and an inequality gi,2{xi,X2,t) 7^ for X2 for 
suitable differential polynomials fi.2, gi,2- We underline that the latter equation and inequality hold on 
Pi. One can treat the system f^2 = 0, (7^2 7^ as the conditions on X2 with the coefficients being some 
differential polynomials in X\ (satisfying Pj). 

Continuing in a similar way, the algorithm produces a triangular system of differential polynomials 
/i,3(a;i,a;2,a!3,*), 9i,a(%i,X2,xa,t),...,fi tn (xi,...,x n ,t), g i<n (xi, . . . ,x n ,t). Thus, at the end x n satisfies 
(on Pi) the equation fi t n(xi, . . . , x n , t) =0 and the inequality gi, n {x\, ■ ■ ■ , x n , t) 7^ treated as a system 
with the coefficients being differential polynomials in xi, . . . , x n -i. 

In other words, suppose that one has a device being able to solve an ordinary differential system 
f(x) = 0, g(x) 7^ in a single indeterminate x. Then the algorithm would allow one to solve the system 
(22 1 consecutively: first producing x\ satisfying fi t i(xi,t) — 0, gi^{x\,t) 7^ 0, after that producing X2 
satisfying fi^(xi,x 2 ,t) = 0, ^2(^1,^2,*) 7^ and so on. 

This completes the description of the output of the algorithm. Now we turn to the issue of its 
complexity. One can bound the orders of the differential polynomials ord(/j jS ), ord(gi, s ) < r-2" := R, 1 < 
i < M, 1 < s < n, the number of the elements in the partition and the degrees M, deg(fi iS ), deg(gi iS ) < 
{Nd f := Q. Finally, the bit -size of the integer coefficients of fi_ s , gi. s and the complexity of the algorithm 
can be bounded by a certain polynomial in I, Q. 

Thus, the number n of the indcterminates brings the main contribution into the complexity bound, 
which is triple exponential in n. Of course, the above bounds have an a priori nature: they take into the 
account all the conceivable possibilities in the worst case, whereas in practical computations considerable 
simplifications are usually expected. 

This illustrates the gain that one can obtain from using tropical geometry to guide model reduction 
and obtain systems with smaller numbers of indeterminates. 



9 Conclusion. 

Tropical geometry offers a natural framework to study biochemical networks with multiple timescales and 
rational reaction rate functions. First, and probably most importantly, tropicalization can guide model 
reduction of ODE systems. We have shown that the existence of quasi-equilibrium reactions and of quasi- 
stationary species implies the existence of sliding modes along the tropical variety. Conversely, when the 
tropicalization has sliding modes along hyperplanes defined by the equality of two monomials, we propose 
an algorithm to decide whether the system has quasi-equilibrium reactions or quasi-equilibrium species. 
This distinction allows correct model reduction, and represents an improvement of methods proposed in 
|RCZL08| . 

The tropicalization represents an abstraction of the ODE model. This abstraction may be not sound 
for some dynamic properties, but may conserve others. If the trajectories of the ODE model are either very 
far or very close to the tropical varieties, they tend to remain close to the trajectories of the tropicalization 
for a while. However, the quality of the approximation is not guaranteed at finite distance from the tropical 
variety. For instance, the exit of tropicalized trajectories from a sliding mode tends to be delayed, and 
smooth trajectories leave earlier neighborhoods of tropical varieties. The example studied in this paper 
also illustrates some properties of bifurcations of the tropicalization, that we have tested numerically. The 
tropicalization qualitatively preserves the type and stability of attractors, but can also introduce delays of 
a Hopf bifurcation. Thus, the tropicalization can only roughly indicate the position of the bifurcation of 
the ODE model. Furthermore, for Hopf bifurcations, the amplitude of the oscillations behaves differently 
for the ODE model and for the tropicalization. In fact, Hopf bifurcations are continuous for the ODE 
model and discontinuous for the tropicalization. 

The tropicalization provides in the same time a reduced model and a "skeleton" for the hybrid dynamics 
of the reduced model. This skeleton, specified by the tropical varieties, is robust. As a matter of fact, 
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monomials of parameters are generically well separated |GR08j . This implies that tropicalized and smooth 
trajectories are not that far one from another. Furthermore, because the tropicalized dynamics is robust, it 
follows that the system can tolerate large relative changes of the parameters without strong modifications 
of its dynamics. 

The dynamics of the model studied in this paper is relatively simple: it has a limit cycle embedded in 
a two dimensional invariant manifold. As future work we intend to extend the approach to more complex 
attractors, such as cycles in dimension larger than two and chaotic attractors. Methods to compute 
tropical varieties in any dimension are well developed in tropical algebraic geometry [BJS + 07] . Given the 
tropical variety, the existence of sliding modes can be easily checked and the pruned polynomials defining 
the fast dynamics calculated. This should lead directly to identification of quasi-equilibrium reactions and 
quasi-stationary species, without the need of simulation (replaces Step 1 in the Algorithm 1). Proposing 
simplified descriptions of the dynamics of large and imprecise systems, tropical geometry techniques could 
find a wide range of applications from synthetic biology design to understanding emerging properties of 
complex biochemical networks. 
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